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^ : ABSTRACT 



Using the two-point Edgeworth series up to second order we construct the weakly nonhnear 
conditional probability distribution function for the density field around an overdense region. 
This requires calculating the two-point analogues of the skewness parameter 5*3. We test the 
dependence of the two-point skewness on distance from the peak for scale-free power spectra and 
Gaussian smoothing. The statistical features of such conditional distribution are given as the 
values obtained within linear theory corrected by the terms that arise due to weakly nonlinear 
evolution. The expected density around the peak is found to be always below the linear prediction 
while its rms fluctuation is always larger than in the linear case. We apply these results to the 
spherical model of collapse as developed by Hoffman & Shaham (1985) and find that in general 
the effect of weakly nonlinear interactions is to decrease the scale from which a peak gathers mass 
and therefore also the mass itself. In the case of open universe this results in steepening of the 
final profile of the virialized protoobject. 
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1 Introduction 



The simplest deterministic model of structure formation proposed by Gunn k, Gott (1972), called 
the spherical model, described the evolution of an ovcrdense region in the otherwise unperturbed, 
expanding Universe. It was extended by Gott (1975) and Gunn (1977) to apply to the evolution 
of matter around an already collapsed density perturbation superposed on a homogeneous back- 
ground. The main prediction of the model (called the spherical accretion or the secondary infall 
model) was that the matter collapsing around the perturbation should form a halo with r~^/^ 
density profile. 

It is much more realistic to assume that the progenitors of structure were not the collapsed 
perturbations but the local maxima (rare events) in the density field which had initially Gaussian 
probability distribution. This was the approach of Hoffman & Shaham (1985) (hereafter HS) who 
applied the secondary infall mechanism to the hierarchical clustering model. They assumed that 
the density peak dominates to some extent the surroundings causing the collapse of the material 
that is gravitationally bound to it. The initial density profile around the peak was approximated 
by the mean density predicted by the initial probability distribution which was assumed to be 
Gaussian. Thus a link was established between the statistical nature of the fluctuations and the 
deterministic character of the spherical model. HS considered scale-free initial power spectra of 
fluctuations and found that the flnal profiles of halos depend on the spectral index n and on the 
density parameter J7. 

One of the key assumptions underlying the calculations of HS was that the matter influenced 
by the peak collapses onto it undisturbed by the background. This is equivalent to the statement 
that the peak identified with some resolution (smoothing) scale collapses while the surrounding 
density field is still linear i.e. its rms fluctuation at this scale is much smaller than unity. Since 
the rms fluctuation grows with time and decreases with the smoothing scale, this approximation 
might be true for very early stages of evolution or very large smoothing scales. The example 
of the Virgo supercluster (which has not yet collapsed) shows that this is not always the case: 
assuming the power spectrum P{k) oc A;~^/^ and knowing that at present the rms fluctuation is 
of the order of unity at the smoothing scale of 10 Mpc, the size of the supercluster being around 
30 Mpc, one can easily estimate that the rms fluctuation at the scale of the supercluster is well 
in the weakly nonlinear regime. 

The purpose of this paper is to present a generalization of the calculations of HS to the case 
of density peaks collapsing in the weakly nonlinear background. In this way we hope to account 
properly for the weakly nonlinear transition between the linear and strongly nonlinear phase of the 
evolution of the perturbation which lacked in the approach of HS. Such a generalization involves 
constructing the weakly nonlinear probability distribution function (PDF) of density around an 
overdense region. The properties of the one-point weakly nonlinear PDFs were discussed by 
many authors (e.g. Bcrnardeau & Kofman 1995; Juszkiewicz et al. 1995) and all studies confirm 
that the weakly nonlinear PDFs develop features which are absent in the linear phase of the 
evolution e.g. the skewness. Those functions describe the density field at a randomly chosen 
point. Here we impose the condition that the point be chosen in the vicinity of a significantly 
overdense region which requires constructing first the two-point weakly nonlinear PDF. Then the 
PDF in this point is given by the conditional probability obtained from the two point PDF with 
the restriction that the density in the second point (the location of the peak) is known and equal 
to a constant. The mean density obtained from this weakly nonlinear PDF is then taken as the 
initial condition for spherical collapse. 
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The reliability of this approach rests on the assumption that the influence of the neighbouring 
fluctuations can be restricted to the weakly nonlinear phase with its only outcome in the form of 
a changed 'initial' density profile that can then evolve independently of surroundings, according 
to the spherical model. Although it is quite obvious that the evolution of a strongly nonlinear 
object is governed mainly by its own gravity and is little affected by its surroundings, it is very 
difficult to determine precisely the moment when such a situation takes place, i.e. a moment 
when we should pass from the weakly nonlinear statistical description to the strongly nonlinear 
deterministic spherical model. 

The paper is organized as follows. In Section 2 we construct the bivariate Edgeworth series 
and calculate the two-point skewness parameter for scale-free power-law spectra and Gaussian 
smoothing. In Section 3 we calculate the conditional probability distribution of the density field 
around a peak and discuss its properties. Section 4 provides the application of the results to the 
spherical collapse model. The concluding remarks follow in Section 5. 



2 The bivariate Edgeworth series 

We consider density contrast field which initially has Gaussian distribution with zero mean. The 
field measured at a given point will be denoted by the symbol 5 while the one measured at the 
distance r from the first point will be called 7. They can be treated as two variables, two fields 
in the same space. Since the fields can in general be smoothed with filters of different scales their 
variances can be different: {S"^) = cj^ and (7^) = r^. The (auto)correlation function of these two 
fields is given by 

(5(x)7(x + r))=e(r). (1) 

Hereafter we will consider the normalized density contrast fields ^ = S/a and v = 7/t. Their 
variances are now equal to unity and their correlation is now given by g = ^joT. The quantity 
Q will be referred to as the correlation coefficient. The joint probability distribution of the two 
variables in the Gaussian case is 

P(/^> ^) = 7, — /! o /(/^' Q) (2) 
27rVl - Q 



where 



f{fi,u, q) = exp 



2(1 - e 



2^ 



(3) 



If the fields were uncorrelated {g = 0), p{fj,, v) would be just a product of two Gaussian distribu- 
tions of /i and V. 

The purpose of this section is to generalize the joint probability distribution function in the 
Gaussian case to the case when the density fields are weakly nonlinear and therefore departing 
from Gaussianity. If the rms fiuctuations of the fields are small (below unity) the fields can be 
expanded around their linear values 5i and 71 respectively 

5 = 6i + 52 + 6'i + --- (4) 
7 = 71 + 72 + 73 H • (5) 

As discussed by Juszkiewicz et al. (1995) the so-called Edgeworth series provides a good 
approximation to the one-point PDF, except for the very tails of the distribution. Up to the 
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third order approximation the Edgeworth series reads 



-1.2/2 



(6) 



where 6*3 and S4 are respectively the third and fourth normahzed cumulants of the density 
contrast field, the skewness and the kurtosis (see e.g. Bernardeau 1994; Lokas et al. 1995) and 
Hn{y) is the Hermite polynomial of the n-th order generated by 



(-1)" 



d" 



di/" 



(7) 



In the construction of the bivariate Edgeworth series we follow the work of Longuet-Higgins 
(1963, 1964) who considered modified Gaussian distributions of weakly nonlinear variables and 
applied them to statistical theory of sea waves. In the second order approximation we have 



1 + 7:('^30-f^30 + 3A2li?21 + 3Ai2-f/^12 + Ao3-f^3o) 
D 



(8) 



where f{fi,h',Q) is given by equation The correlation coefficient g, the Hermite polynomials 
Hmn and the two-point moments Xmn will be discussed in the following. 



2.1 The correlation coefficient 

The autocorrelation function of density field measured at two points separated by distance |r| = r 
can be calculated using the relation 

Ur) = J j^P{k)W\kR) e^*^- (9) 

where P{h) is the power spectrum of the density fluctuations and the field is smoothed with a 
window function W of radius R which in what follows will be assumed to have a Gaussian form 

W{kR) = e-'='«V2_ (10) 

In what follows we will restrict the calculations to the case when the field at both locations has 
the same properties and is smoothed with the same smoothing radius so that the variances are 
equal: a = t. 

For the scale-free power spectra 

P{k) = Ck"", -3 < n < 1 (11) 

we obtain 

M:) _ r(i + i) ^1/2 (l) ^_../4 .^2) 

^- ^2 - 2 r(^) "/2l^4 ) ' ^^^^ 

where c = r/R and L'^{x) are Laguerre polynomials. Here we have used the fact that for the 
scale-free power spectra ( [lT|) and Gaussian smoothing (|lO| ) the linear variance of the density PDF 
is given by 

= i^i) = D\t) j j^P{k)W\kR) = CD\t) ^^^yJ^^ . (13) 
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Hereafter the variance (5^) will be always approximated by the linear value (^) which is sufficient 
in the construction of the two-point Edgeworth series up to second order. The result (12) takes 
simpler form in terms of the degenerate hypergeometric function 

Q{n,c) = iFi I 1 . (14) 

We see that the correlation coefficient depends on the scales r and R only via their ratio c. Upper 
panel of Figure || shows the shape of g{c) for different power spectrum indices n = —2, —1.5, — 1. 

2.2 The bivariate Hermite polynomials 

The quantities Hmn in (§) ai'e two-dimensional analogues of the Hermite polynomials Hn appear- 
ing in the one-point Edgeworth series (^) that were generated by equation The bivariate 
Hermite polynomials can be calculated using similar formula 

^, e) = /(m, ^, g)H^nif^, u, q). (15) 

The polynomials of orders needed in this calculation can be written in the form 

1 

{l-Q 



where 



(m - QJ^f of ^ 
'ISO = — i 2 ^l/" ~ 

h,, = 2,(;. - H - - ^?/^) + ^^^^T^^^^f^ (17) 



hi2 = 2g{u - gfi) - (/i - gj^) + 

{v- guf 

ho3 = —. n 2>(u-gfi). 

I — g^ 



{v - gfi){fi - gvf 
1-^2 

(// - gv){v - g^if 



,2 



2.3 The two-point skewness 

The coefficients \ij in (P) are defined by the reduced moments (cumulants) Kij in the following 

way 

A,, = (18) 

while the cumulants of the fields 5 and 7 are given by the connected part of the moments 

Kij = {d'^')c- (19) 

Using the perturbative expansion of the fields, equations (@)-(|5|), we find up to second order 

K30 = {Sf) + 3{5l52) 

K21 = {Sh2) + 2{5i52ii) (20) 
K.12 = (^27i) + 2((5i7i72) 
K03 = (7i) + 3(7^72). 
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Figure 1: The values of the correlation coefficient q, the two-point skewness iS'12 and the ratio 
S12/Q as functions of the distance between the two points (in units of the smoothing scale), 
c = r/R, for different power spectra. 
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The moments K20 and K02 are just the variances (5^) and (7^) respectively, which, as stated 
before, can be both approximated by the Unear value o"^ given by equation (|T^). 



The moments K30 and kqs are related to the ordinary one-point skewness parameters (see e.g. 



Peebles 1980; Juszkiewicz et al. 1993). For Gaussian fields (df) = 0, and we have 

K30 = 3{Sj52) = Ssa^ 



(21) 



where 5*3 is the dimensionless skewness parameter which for a Gaussian filter depends on the 
spectral index n in the following way (Lokas et al. 1995) 



3 2F1 



n-|-3n-|-3 3 1 
2 ' 2 '2' 4 



8 

n+- 



2F1 



n-|-3n-|-3 5 1 
2 ' 2 '2' 4 



(22) 



The other one-point moment kqs = S^t'^ but since we have already put r = o" we have K03 = 1^30 ■ 

The two-point moments K21 and K12 can be calculated using the same method as in the case of 
one-point skewness (Lokas et al. 1995) but now we have to keep the dependence of the moments 
on distance r = |r|. The two-point moments also scale as cj^ so it is convenient to calculate them 
in the normalized form. Using equation (13) we have for K21 /cr'^ 



1 



4487r6i?2{n+3)^4 

P(p)P(Q)J(p + q,p,q; 



d?p / d^q l^(p)VF(g)VF(|p + q|)e^(P+^) '" 



(23) 



and 



2(<5i527i>(0 



2247r6i?2(n+3)^4 

P(p)P(g)J(p + q,p,q; 



d?p / d?q Ty(p)VF(g)'W^(|p + q|)e'P-'" x 



(24) 



where J is the kernel of the second order perturbative solution for density field. The integrals 
(23) and ( |2^ ) will be evaluated for scale- free power spectra (pT|) and Gaussian smoothing ([lO|). 
The simplest way to perform the integrals in ( [2^ ) is to change variables from p, q to 1, p so that 

p + q 



1 and 1 • p = / p/3. Together with the expression for the variance (13) this produces 
{^h2){c) 2 



7r2( 



n+3' 



dp j dl sin(cO I e" 



X / d/3 J(/, p, [5) {f + p2 _ 2Zp/3)"/2 



(25) 



where c = r / R as before. The integration over the angular variable /? can still be done analytically 
for integer or half-integer values of the spectral index n. The remaining integrations are done 
numerically. In the case of (|23) the integration over angular variables gives 



2(5i527i)(c) 



8^^ 



^r2(n±3) c 

^^1/2 (W) - 



Jdpjdq sm{cp) qn+3/2 



P Q 

- + -] h/iipi) + ^4/2 (P?) 



(26) 



By using the expansion of Bessel functions 

00 



1 



m=0 



m\r{u + m + l) \2 



u+2in 



(27) 



7 



and the facts that 



f°° 2 

/ q'^e-'^ dq 
Jo 

I p e ^ sin(cp) dp 
Jo 



ir 

2 



a + 1 



- r 1 + - 
2 V 2 




for b> -2 



(28) 
(29) 



we obtain the result in the form of a series of combination of gamma and degenerate hypergeo- 
metric functions or Laguerre polynomials 



2{6i627i){c) 



2r2(n±3) 



-eV4 g r(^+m) ^ 
m=0 



! 22^^ 



34r(n±2+m) ^1/2 



21 r(m + |: 



n/2+m 



2r(m+|) 



i I h m 



L 



1/2 

n/2+m+l 



+ 



+ 



+ TTi i h m 



2 / V 2 

(B±3 + m)r(^ + m) 



L 



1/2 

n/2+m 



21 



4r(m + 



L 



1/2 

n/2+m+l 




(30) 



The values of 



'21 



K21 
^4 



(31) 



will be hereafter called the two-point skewness parameter. Because of the symmetry between the 
two locations we obviously have = 5*21. In the following calculations we will therefore use 
only the symbol Si2. The numerical values of S12 for different values of c and n = —2, —1.5, — 1 
are given in Table |l| and plotted in the middle panel of Figure [l|. 

The lower panel of Figure |^ shows the ratio S12/ Q- As discussed by Bernardeau (1996) for the 
case of top-hat filter such quantity should approach a constant at c ^ 00. Figure |l] shows that 
the same behaviour is observed for Gaussian filtering. Using diagrammatic representations of the 
two-point moments Bernardeau (1996) has shown that the part of 5i2 given by equation ( p3|) 
behaves like in this limit and therefore is negligible compared to the part given by equation 
(24) which behaves like q. For large c the value of the sum in (|30| ) is perfectly approximated by 
the m = element of the series and we find 



5*12 = 512^1 



(32) 



where 



(33) 



47 n 
- 21 " 3 • 

It is worth noting that the moment S12 given by equation ( |33| ) is exactly equal to the corresponding 
quantity C2,i calculated by Bernardeau (1996) for scale-free power spectra and top-hat smoothing. 



For n = —2, —1.5, —1 the numerical values of si2{n) are respectively: 2.90, 2.74 and 2.57. It is 
seen in Figure |l] that the values of S12/Q indeed approach them for large c. The value of 5i2 is well 
approximated by the one given by equation (^) within 10% for c > 4 (and the approximation 
works much better for n = —1 and n = —1.5 than for n = —2). In the following, however, we will 
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5'l2 

c = r/R n = —2 n = —1.5 n = — 1 






4.022 


3.714 


3.468 


0.01 


4.022 


3.714 


3.468 


0.1 


4.018 


3.708 


3.460 


0.5 


3.913 


3.563 


3.282 


1 


3.614 


3.159 


2.790 


2 


2.742 


2.062 


1.547 


3 


1.957 


1.212 


0.7212 


4 


1.441 


0.7545 


0.3659 


5 


1.126 


0.5226 


0.2206 


6 


0.9237 


0.3914 


0.1494 


7 


0.7830 


0.3081 


0.1084 


8 


0.6796 


0.2510 


0.08233 


9 


0.6003 


0.2097 


0.06471 


10 


0.5377 


0.1787 


0.05222 



Table 1: The values of the two-point skewness parameter S12 for density field as a function of 
the spectral index n and c = r/R. 



be also interested in the region of c of the order of few, where this approximation is not sufficient 
and cannot be adopted. 

In the limit of c — > 0, which is equivalent to r ^ 0, the two points at which we measure the 
density field converge and we expect that the two-point skewness reproduces the one-point value 
of 53 given by equation As Table 1 and Figure |l| prove this is indeed the case. In this limit 
the two parts (|2^ ) and ( P4^ have comparable contributions: S-^/?) and 253/3 respectively. 

Using the results of this section the bivariate Edgeworth series (^) can be rewritten in the 
following way 

P{l^,^) = 7;— 7==f ^) (34) 



27rv/T - Q 



|l + g^j^ - g2)2 ['^'3(^30 + ^03) + 3S'l2(/ll2 + ^2l)]| 



3 Conditional probability distribution around a peak 

Let us now suppose that the value of one variable is known: v = a. The conditional probability 
distribution of the other variable is defined as 

1/ = a) = — — (35) 

J_oo P(M,a) An 

where the nominator is the bivariate distribution with the second variable put equal to a constant 
and the denominator is the marginal distribution of the variable u = a. 

In the case of two-point Gaussian distribution (pt) the resulting conditional distribution for 
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reads 



p{n,u = a) 



1 



exp 



(/i - gg)^ 
2(1 - g^) 



(36) 



V2vr(l - 

This is again a Gaussian but the variance was changed from unity to 1 — and the mean value 
of fi was moved from zero to fi = ga. 



In the case of two-point Edgeworth series ( p4[ ) the conditional distribution of //, provided the 
value of u is known and equal to o, is obtained from equation (|3^) using the marginal distribution 
of v 

p{fj,, u = a) d^i = 



-aV2 



27r 



6 



3o) 



(37) 



which is just the one-point Edgeworth series (|6|) up to second order with u = a. Since only the 
lowest order correction in the Edgeworth series was introduced the result must be expanded in 
powers of a and only the term linear in a should be kept 



1 + 



a 



exp 



in - gaf 



X 



2(1 - g-^ 

[-5*3(^30 + /lOs) + 3S'l2(/ll2 + /l2l)] 



(38) 



6 



3a) 



where hr, 



6(1 - g^ 

hmni^J'l g) are given by equations (|D 

First let us check the behaviour of the conditional probability ( |38|) in the limit of the corre- 
lation coefficient g ^ i.e. when the two considered points are separated by large distance and 
therefore uncorrelated. According to the results of the previous section in this limit also S12 — > 
and the distribution (H) becomes 

1 



p{fi,u = a)|g_o 



V/2 



3^) 



(39) 



which is just the one-point Edgeworth series (P), as expected. The probability distribution for 
the density field is not supposed to depend on a particular value of density contrast at infinite 
distance. 



In order to test the behaviour of the distribution ( pq ) for different parameters g, a and 
a we calculate its lowest order moments with respect to the mean value e (the subscript 
E will identify the quantities obtained for the distribution (|3^)) and compare them with the 
corresponding quantities of the Gaussian conditional distribution (p6D. Let us recall that for 
Gaussian distribution we found (in agreement with Dekel 1981, Peebles 1984 and HS) 



((/^-(^))')g 



g a 
1 - 
0. 



After tedious but straightforward calculations for the Edgeworth conditional probability 
get 



a 



g a + -(a 



1){Si2-qS^) 
g^ + a a {S12 - 2gSi2 + g'^S^ 



a {S3 - ^gSi2 + 3g'^S_ 



(40) 
(41) 
(42) 

I) we 



(43) 

(44) 
(45) 



10 



From equation ( |4^ ) we also have the dispersion to the lowest order in a 

Clearly the characteristics of the conditional Edgeworth distribution are given by the Gaussian 
ones plus correction terms proportional to a and a function of a (except for the third moment 



which is independent of a) . As shown earlier in equation ( 3£ ) in the limit of — > the conditional 
Edgeworth distribution approaches one-point Edgeworth series. This must also apply to its 
moments. Using the fact that at small q the values of two-point skewness S12 vanish we again find 
that in this limit the third order moment ( p5[ ) approaches the value of the one-point Edgeworth 
series, S'^o. The same is true for the lower moments (|4^ ) and (^4|). 



An independent check of the results (|4^)-(|4q) is provided in the limit oi q ^ \. Then the 
distribution p(/i,z^ = a) should approach the Dirac's delta function, boi^ii — a), with the mean 
equal to a and the second and third moment equal to zero. Recalling that in this limit 512 S3 
such results can be immediately reproduced from equations (^3|)-(^). 



An interesting application of these results from the point of view of the theory of structure 
formation is to study the density distribution around an overdense region. Although the results 
presented so far apply to arbitrary values of a we will from now on focus on regions which are 
able to dominate their surroundings, that is their density contrast is bigger than one standard 
deviation. This corresponds to assuming a > 1. A region of such overdensity will be called a 
peak although it might not be a maximum in a strict mathematical sense. However, we expect 
that points of a ^ 1 most probably correspond to local maxima. As reasoned by HS, this is also 
true for mild values of a ^ 1 (see also Adler 1981). 



If a > 1 the moments (|40|)-(|45D can be interpreted as characteristics of the density distribution 
around a density peak. The Gaussian values are the well known results of linear theory while 
the values for the Edgeworth conditional distribution provide the corrections introduced by the 
fact that the field on a given smoothing scale has become weakly nonlinear. The direction of the 
effect depends on the numerical values of the one-point and two-point skewness parameters and 
the correlation coefficient. 

It is clear from Figure |l| that for the spectral indices considered here {n = —2, —1.5, —1) and 
nonzero distance from the peak we have S12/ Q < S3 which can be rewritten as 

S12 - 0S3 < 0. (47) 



This proves that according to equation ( p^ ) for a > 1 the correction to the mean normalized 
density with respect to the Gaussian value is always negative. Figure |2| shows the expected nor- 
malized density contrast (/i) for the Edgeworth (|4^) and Gaussian ( [40| ) conditional distributions 
as a function of the distance from the peak (in units of the smoothing scale), c = r/R. Each 
panel shows results for different scale-free power spectrum with spectral index n = —2, —1.5, —1. 
In each panel the Edgeworth results are plotted as thicker solid lines while the thinner dashed 
lines show the Gaussian ones. The four pairs of lines in each panel correspond from bottom to 
top to the four chosen values of the peak's height: a = 1, 2, 3, 4 (for a = 1 the results are the same 
in both Gaussian and Edgeworth cases). The Edgeworth results are always below the Gaussian 
values so it is clear that the effect of weakly nonlinear interactions is to decrease the expected 



density around an overdense region. Equation (43) shows explicitly that the effect grows with a 
and fj. 
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Figure 2: The mean normalized density contrast as a function of c = r/R. The dashed hues 
correspond to the Gaussian case while the solid ones show the results obtained from the Edgeworth 
approximation with a = 0.2. Each panel shows results for a = 1,2,3,4 and different power 
spectrum. 
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Juszkiewicz et al. (1995) have tested the third order Edgeworth approximation against N- 
body simulations and found that it is accurate for density contrast up to (5 = acr = 1. Here we 
are using only the second order approximation so the range of validity of these results is probably 
even more restricted. Therefore in the following it will be assumed that aa is always below unity. 
All the Edgeworth results in Figure |2| are given for a = 0.2. At this value of the rms density 
fluctuation we still have aa < 1 even for o = 4 so we can expect the Edgeworth series to be a 
good approximation. 

In the case of the variance the weakly nonlinear corrections work in the opposite direction: 
their effect is to increase the value of the variance (or dispersion) with respect to the linear case, 
because for the power spectra considered here (and nonzero distances from the peak) we always 
have 

Si2 - 2gSi2 + g^Ss = (1 - g)Si2 + gigSs - 812) > (48) 



where the inequality (47) was used. The effect grows linearly with a and a as equation (| 



states. In Figure ^ we plot the dispersion of the conditional Edgeworth distribution (^) and the 
Gaussian dispersion which is independent of a and equal to — The quantities are shown 
as functions of c = r/R for a = 0.2, a = 1, 2, 3, 4 and different power spectra. As in Figure Q the 
thicker solid lines in each panel correspond to the Edgeworth results and the thinner dashed lines 
to the Gaussian ones. Although the Gaussian dispersion is always below unity and approaches 
it at large c, the Edgeworth values reach a maximum, which is above unity at c of the order of a 
few and then decrease down to unity at large distances. 

An important conclusion coming from the behaviour of weakly nonlinear expected density 
and its variance is the following. It is clear that a local density peak that rises significantly above 
the noise should gravitationally dominate its surroundings out to some distance. A reasonable 
measure of the distance, up to which a coherent structure around the peak is expected, is the 
scale r^oh at which 

(A.) = ((/i-(M))2)^/2. (49) 

Because of smoothing it is more convenient to measure the distance r in units of the smoothing 
scale R so in what follows we will use Ccoh = fcoh/R instead of rcoh- The mean (/i) as a function 
of c up to Ccoh may be therefore treated as the density profile of matter that is bound to the peak 
in a statistical sense. 

In the Gaussian case equation (^9|) leads to g{n,Ccoh) = 1/ Va^ + 1 which can be solved 
numerically for Ccoh given the shape of g, equation (0). Finding Ccoh in the case of Edgeworth 



approximation requires fitting the function 512(0) and then equation ([49[) can also be solved 
numerically. The results for both cases are shown in Figure ^ as dependent on the height of the 
peak o, for three different values of cr = 0.05, 0.1, 0.2 and three spectral indices n = —2, —1.5, —1. 
Since the mean density is decreased and the dispersion is increased the effect of nonlinearity 
is always to decrease the coherence length. The magnitude of the effect grows with the rms 
fluctuation a. Therefore the density field around a peak embedded in the weakly nonlinear field 
decorrelates much faster with distance than in the linear case. 

Note that the scale of each panel in Figure ^ is different: the values of Ccoh are generally 
significantly larger for lower spectral indices but this reflects mainly the amount of large-scale 
power in the initial linear power spectrum. To take this into account the weakly nonlinear 
coherence length should be compared to the corresponding Gaussian one. Then the dependence 
on the power spectrum is clear: for a given a and a the change in the value of Ccoh with respect 
to the Gaussian one is larger for lower spectral indices. For example for a = 3 and a = 0.1 the 
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Figure 3: The dispersion of the normahzed density contrast as a function of the distance from the 
peak, c = r/R. The dashed hne in each panel shows Gaussian results which are independent of the 
peak's height while the four solid lines correspond to the results obtained from the Edgeworth 
approximation with a = 0.2 and a = 1,2,3,4. Each panel shows results for different power 
spectrum. 



14 




Figure 4: The length of coherence (in units of the smoothing scale) defined by equation (|49| ) as a 
function of the peak's height, a. The dashed line in each panel shows the results in the Gaussian 
case. The solid lines correspond to the results obtained from the Edgeworth approximation with 
a = 0.05,0.1,0.2, with larger a producing a curve that departs more from the Gaussian one. 
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weakly nonlinear Ccoh is decreased with respect to the linear value by 11%, 14% and 21% for 
n = —1, —1.5, —2 respectively. 

In general the coherence length grows with the height of the peak, a, which is reasonable: 
higher peaks should dominate their surroundings to farther extent. However, for largest values 
of (T = 0.2 we observe that the curves in Figure ^ have a maximum at o ~ 5. Since in this case 
acj = 1 it is probably the effect of breakdown of the Edgeworth approximation. 



Another important quantity to characterize the conditional distribution ( |38D is the third 
moment (^). For the Gaussian distribution this moment is always equal to zero, while for the 
weakly nonlinear case a good measure of it (and the asymmetry of the distribution) is the quantity 

{{ll - {ll)f)E CJ (^3 - 3^>5i2 + 3i>25i2 - Q^S; 



(50) 

where we have used equations (Q)-(^5p, the quantity was expanded in powers of a and only the 
term linear in a was kept. 

Figure |5| shows the values of the normalized third moment ( [50| ) as a function of c for cr = 0.1 
and different power spectra. It is obvious that the shape of the curves would be the same for 
any other a. It is worth noting that the value (50) which provides the lowest order correction to 



the Gaussian case does not depend on the height of the peak, a. The third moment normalized 
in such way grows with c for all considered spectral indices. For large c the values (at a given a 
and c) are higher for lower spectral indices and they approach those of the one-point Edgeworth 
series, S30", with 53 given by equation (p^: 4.022, 3.714 and 3.468 respectively for n = —2, —1.5 
and —1. It is interesting that up to c ~ 3 the dependence on the spectral index is weak but 
opposite: the third moment grows with n. 



4 Application to spherical collapse 

The dynamical evolution of matter at the distance Cj from the peak is determined by the mean 
cumulative density perturbation within Cj which is given by 

Ai(Q) = 4 r '5(c)c2dc (51) 



cf Jo 

where we put 5(c) = (j(^) with (^) given by equations (l40|) and (^3|) in Gaussian and Edgeworth 
case respectively. In the linear case the function Aj(cj) is well approximated by 

\ ao for Cj ^ 1 

Ai,G(Ci) = \ ... -(n+3) . ^ T (52) 

y aan[n)c^ for Cj 1 

where h{n) is a numerical factor. Its values for different spectral indices are given in Table ^. 
The approximation at large q is accurate to within 10% for Cj > 5. 



In the weakly nonlinear case equations (|4^) and (^) give 

j acj for Cj ^ 1 

^'^^''^^ = 1 aah{n) [l - a{S3 - si2)ia^ - l)/(2a)] cr^"^^) fo^ ^ (53) 



where ^3 and S12 are given by equations ( |22D and ( pSD respectively. The factor 5*3 — S12 is of 
the order of unity; the exact numerical values for different n are given in Table ^. Equation (53) 
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Figure 5: The properly normalized third moment of the conditional distribution as a function of 
the distance from the peak, c = r /R, for a = 0.1. Different lines show results for different power 
spectra. At large distances the curves flatten to reach the limiting values of S^a. 
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n 


h 


S3 - S12 


-3.0 


1. 


1.619 


-2.5 


1.659 


1.329 


-2.0 


2.659 


1.117 


— i.o 




U.y 10 1 


-1.0 


6. 


0.8966 


-0.5 


8.296 


0.8747 


0.0 


10.63 


0.9064 


0.5 


12.27 


0.9897 


1.0 


12. 


1.124 



Table 2: The values of the parameters h and S3 — S12 as functions of the spectral index n. 



shows that in the limit of large Cj weakly nonlinear interactions do not change the slope of the 
average initial density profile. Instead the effective height of the peak is changed: the profile ( |53[ ) 
at large Cj can be viewed as a linear profile (^) with the effective height of the peak 



^eff 



1 



2a 



(54) 



If we now assume that the matter of average density Aj within distance from the peak 
collapses undisturbed onto the peak, the spherical model can be applied. In what follows we 
compare the predictions of the spherical model with Gaussian initial conditions to the ones 
obtained with the initial conditions settled by the Edgeworth approximation. 

According to the spherical model the radius of maximum expansion is related to the initial 
radius of the shell rj in the following way 



A, + l 



(55) 



where 6c = — 1 and Qi is the density parameter at some initial epoch. The maximal radius 
is related to the radius of the shell after virialization, r, in the following way 

r = frm, (56) 

where / is of the order of 1/2. 

From now on we will focus on the analytically tractable case of large Cj where the correspond- 



ing Aj are much less than unity. Combining (55) and (56), expressing distances in units of the 
smoothing scale R, r = cR, and expanding in powers of Aj we have 



1 A. 



1 



/ A, + l yMc + (l + '5c)A,, + 0(A2)]. 



(57) 



The expansion (p^) is in fact an expansion in powers of c. 



-(n+3) 



small parameter but it is 



equivalent to the expansion in a if the terms of the order of Cj ^^"'"'"^^ and higher are neglected. 

The value of 5c defines the radius cq of the shell, the total energy of which vanishes. All shells 
of radii Cj < cq are gravitationally bound to the peak and will eventually collapse onto it. The 
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condition of vanishing energy leads to Aj(co) = 5c and in the hnear case we have 



Co,G 



while in the weakly nonlinear approximation 



C0,E 



CO,G 



cTjSs - su){a^ - 1) 
2a{n + 3) 



(58) 



(59) 



The weakly nonlinear cq is therefore decreased with respect to the linear one. 

Since all shells which are gravitationally bound to the peak have Ai(ci) > Sc the term (5cAi 
in equation (|5^) can be neglected as being of the order of A? and we have 



c ~ f 



n+3 



(60) 



The density run of the collapsed object can be estimated by assuming that the material 
originally at the shell of radius Cj ends up at c so that 



p{c)c'^dc = pi{ci)cf da. 



The initial density of the shell of radius Cj 



PiiCi) = Pb,i[^ + '^i(cj)]> 



(61) 



(62) 



can be approximated by pb^i, the background density at the initial epoch, since at large Cj con- 
sidered here the expected value of 6i is very small in linear case and in fact even smaller in the 
weakly nonlinear approximation as was demonstrated in equation (^3[). The density profile is 
therefore in general given by 



P{c) 



Pb,i T 



[(co/c, 



\n+3 



(63) 



/; (n + 4)(co/ci)"+3-l- 
In the two limiting cases of interest that is of q much smaller than and comparable to cq we have 



N 3(n+3)/(n+4) 

) for Ci <C Co 

c / 



/ N Pb,i f f <^0Y f ^ 

P{C) = -^T- — for Ci < Cq. 

n + 60c V c / 



(64) 
(65) 



The dependence of p{c) on the ratio cq/c and the form of the two limiting cases are then the 
same in the weakly nonlinear as in the linear case discussed by HS. When co is very large or 
very small the weakly nonlinear corrections will not affect much the linear results and the slopes 
( |6^ and (^) will be preserved. The case of cq — > oo corresponds to = 1. If at present we 
have = 1 then surely in the past also = 1 and all profiles should have the asymptotic form 



(64). If the present universe is open the profile depends on the collapse time of the structure: 
the structures that collapsed later (with lower Oj) should have smaller cq and steeper profiles 
because, as shown by HS by fitting numerically some c~" profiles to the general formula (^) for 
given values of cq, the smaller cq, the steeper is the density profile. 
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n 




Co,G 


CO,E 


ac 


aE 


-2 


0.7 


6.2 


2.7 


2.2 


2.5 




0.8 


10.6 


4.7 


1.9 


2.1 




n n 
U.y 


zo.y 


iU.o 


i. ( 


1 Q 
l.O 




1.0 


oo 


oo 


1.5 


1.5 


-1 


0.7 


3.7 


2.9 


2.8 


2.9 




0.8 


4.9 


3.8 


2.5 


2.6 




0.9 


7.3 


5.7 


2.2 


2.3 




1.0 


oo 


oo 


2.0 


2.0 



Table 3: The slopes of profiles c " fitted in the range 1 < c < 10 in the linear (G) and weakly 
nonlinear (E) case for spectral indices n = —2 and n = —1 and different values of Oj. 



The impact of weak nonlinearity can therefore be seen only at intermediate values of cq when 
the correction to cq given by equation (^) is significant. To estimate the correction to the slope 
introduced by the change in cq we perform fits similar to those of HS here first with cq given by 
the linear value (58) and then the weakly nonlinear one, (^9|). The fits should be treated only as 



indicative because they are not rigorous perturbatively. We adopt the maximum value of acr = 1 
allowed by the Edgeworth approximation and find that the weakly nonlinear cq^e (for reasonable 
n > —2) can be as low as co,g/2. The results of the numerical fits of the form c~" to the formula 
( |63|) in the range 1 < c < 10 for the values of cq corresponding to different values of < 1 in 
the case of n = — 1 and n = —2 are given in Table ^ 

The results clearly display the dependences on Jlj and n discussed by HS: profiles are steeper 
for lower ilj and higher n. The weakly nonlinear approximation however predicts steeper slopes 
for i7j < 1 than the linear approximation does. Only the case of fij = 1 remains unaffected. 
Therefore in general weakly nonlinear corrections act in the same directions as decreasing Qi or 
increasing n. This must be taken into account were the slopes of the profiles used to determine 
these cosmological parameters. 

Another important quantity characterizing the protoobject is its mass. We assume that the 
mass contained initially within radius cq 

4:7V 

M(co) = —{coRfpb,,{l + Ai(co)) (66) 

does not change as the shells collapse and may be treated as the final mass of the protoobject. In 
this case the quantitative prediction concerning the correction to the mass gathered by the peak 
can be made. In the linear approximation the mass Mq is calculated with co,g and Aj^c(co) given 
by equations (|5^ and (|5^) respectively while in the weakly nonlinear case [M^) they should be 
replaced by the values given in equations (^) and (|53|). For the ratio of the masses we obtain 



Mg 2(n + 3)a ^ ^ 

Thus the mass gathered by the peak predicted by weakly nonlinear approximation is always 
smaller than in the linear case and the correction is larger for lower spectral indices. 
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5 Concluding remarks 



The picture emerging from the presented results is the following. In the weakly nonlinear field 
the average overdensity around the peak, which drives the evolution of matter around it, behaves 
as if produced by the peak of reduced height embedded in the linear field. Weakly nonlinear 
interactions decrease the size of region from which the peak gathers mass and the slope of the 
final density profile is steepened. The type of dependence of the profiles on the cosmological 
parameters Q and n derived by HS is preserved. 

It should be remembered that the results presented here were obtained for the mean density 
perturbation, (Aj £;), using the mean normalized density contrast (/i) calculated from the condi- 
tional distribution (|38|). In order to specify the uncertainty of the results it would be desirable 
to know also the standard deviation of Aj.^;, cja, since only below the scale determined by the 
condition Ai^E = the perturbation Ai^E is significant. Such scale would always be finite and 
dependent on the power spectrum. In the case of = 1 universe this would mean that we cannot 
indiscriminately apply the limit cq — > oo, which is rather unrealistic: gravitational influence of 
a peak cannot reach infinite distances because there are always neighbouring peaks that gather 
mass. Calculating cja in the weakly nonlinear regime would however require constructing the 
three-point Edgeworth series, which is beyond the framework of this paper. 

There are other limitations to the weakly nonlinear approach to the collapse of peaks presented 
here. The most obvious one is the condition aa < 1 that must be satisfied for the Edgeworth 
approximation to be valid. Is the evolution in this regime representative of the fully nonlinear 
clustering? Another question is how the linear and weakly nonlinear results should be compared. 
Throughout this paper we assumed the linear value of the variance of the density which is correct 
as far as we use the lowest order Edgeworth approximation. However, Lokas et al. (1996) have 
shown that the value of a, the typical fluctuation of density, itself changes during weakly nonlinear 
evolution. It remains to be understood how the evolution of a typical fluctuation and evolution 
of matter around a peak are related. 

Observational data concerning the halo density profiles are presently very sparse and derived 
most often from the fact that the observed rotation curves of galaxies are flat, which in the 
Einstein-de Sitter universe leads to the profile that according to formula (|6^ ) corresponds 
to n = — 1. Such spectral index is roughly consistent with what is observed at weakly nonlinear 
scales but we can hardly go beyond such rough consistency checks as far as observations are 
concerned. Therefore although such data are of potential cosmological interest the only way to 
verify the theoretical predictions is to resort to N-body simulations. 

The profiles of dark halos have been measured in N-body simulations by many authors in- 
cluding Dekel, Kowitt & Shaham (1981) Quinn, Salmon & Zurek (1986), Efstathiou et al. (1988) 
and Crone, Evrard & Richstone (1994). The results of such simulations seem to converge to a 
statement that the present structure of halos retains information on the initial conditions and 
displays the dependence on cosmological parameters. Crone et al. (1994) performed simulations 
for power-law spectra and different density parameters f^o and confirmed the overall trend of 
steeper density profiles with increasing n or decreasing as predicted by HS. However, the fitted 
slopes were systematically steeper than those given by HS. As discussed in the previous section 
only in the case of the open universe such discrepancy can be assigned to the weakly nonlinear 
effects. The detailed comparison of the predictions of perturbation theory with the results of 
N-body simulations requires further assumptions (concerning e.g. the epoch of the formation of 
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objects, the most realistic value of aa to be adopted) and is currently under way. 

Recently Sheth &; Jain (1997) proposed a new derivation of the slopes of halo profiles from the 
shape of nonlinear correlation function based on the assumption that the halo-halo correlations 
can be neglected i.e. that for sufficiently small scales the input to the correlation function comes 
mostly from one halo. The slope they obtain, a = 3(n + 4)/(n + 5), can be derived within the 
formalism of HS, and would be obtained instead of the result (|6^), if the initial profile of the 
cumulative density is chosen to be Cj ^"''^^^ instead of Cj '^""'"^^ as in equations (|5^)-(|53|). Sheth 
& Jain (1997) claim that such an initial profile is well motivated by the results of Bardeen et 
al. (1986) who found that as the peak height decreases the initial profile becomes steeper than 
Cj ^"^^^ . In this case the approximation used here (i.e. regions of high enough density, and not 
necessarily peaks, are the progenitors of structure) would be less accurate. 

This objection can only be valid under condition that smaller peaks can collapse by themselves 
to form objects. N-body simulations performed by Katz, Quinn &; Gelb (1993) and van de 
Weygaert & Babul (1994) suggest the opposite: there is no clear correspondence between peaks 
in the initial density field and the actual halos for smaller peaks. For example, Katz et al. (1993) 
find that most of the peaks expected to form cluster-size objects (a > 4) indeed end up in such 
objects while the peaks from which the galaxy-size objects form (a ~ 2) may end up in a galaxy- 
size group as well as merge with a larger object. Even more possibilities for the history of a peak, 
including its breaking into few distinct halos, were found by van de Weygaert & Babul (1994). 
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